/*	TargetedAbundnaceAnalysisSeriveImpl.java
 *
 *	Copyright (c) Manateem Challenge Project, Zaia Group, and contributors 2010
 *
 *	This file is part of MSGToolSuite ('this project').
 *
 *	This project is free software: you can redistribute it and/or modify it
 *	under the terms of the GNU General Public License as published by the Free
 *	Software Foundation, either version 3 of the License, or (at your option)
 *	any later version.
 *
 *	This project is distributed in the hope that it will be useful, but WITHOUT
 *	ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 *	FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
 *	more details.
 *
 *	You should have received a copy of the GNU General Public License along with
 *	this project. If not, see <http://www.gnu.org/licenses/>.
 */
package edu.bu.manateem.msg.analysis;

import java.io.File;
import java.util.Collections;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
import java.util.Set;
import java.util.TreeMap;

import edu.bu.manateem.msg.analysis.util.MzXmlUtils;
import edu.bu.manateem.msg.heatmap.HeatMapBuilder;
import edu.bu.manateem.msg.model.Compound;
import edu.bu.manateem.msg.model.CompoundParent;
import edu.bu.manateem.msg.model.IntensityValue;
import edu.bu.manateem.msg.model.MassScan;

import org.slf4j.Logger;
import org.slf4j.LoggerFactory;
import org.systemsbiology.jrap.MSXMLParser;
import org.systemsbiology.jrap.Scan;

/**
 *	<p>
 *	</p>
 *
 *
 *	@author
 *		Christopher Jacobs
 *
 *	@version
 *		1.0
 *
 *	@changelog
 *		Apr 21, 2010	- Version 1.0.
 */
public class TargetedAbundanceAnalysisImpl
{
	/*
	 * TYPES
	 */

	/*
	 * CONSTANTS
	 */

	// Serialize
	private static final long serialVersionUID = 1L;

	// Logging
	private static final Logger LOG =
		LoggerFactory.getLogger(TargetedAbundanceAnalysisImpl.class);

	/*
	 * FIELDS
	 */

	private static double rtSearchWidth =
		TargetedAbundanceAnalysis.DEFAULT_RT_SEARCHWIDTH.doubleValue();
	private static double rtSumWidth =
		TargetedAbundanceAnalysis.DEFAULT_RT_SUMWIDTH.doubleValue();
	private static double mzSumWidthLeft =
		TargetedAbundanceAnalysis.DEFAULT_MZ_SUMWIDTH.doubleValue();
	private static double mzSumWidthRight =
		TargetedAbundanceAnalysis.DEFAULT_MZ_SUMWIDTH.doubleValue();

	/*
	 * CONSTRUCTORS
	 */

	/*
	 * INITIALIZERS
	 */

	/*
	 * METHODS
	 */

	/**
	 * Combines (sums) search windows belonging to a parent compound,
	 * since a compound may have more than one 'child.' 
	 * 
	 * @param parentCompound
	 * @return a mapping of retention time to intensity value
	 */
	public static Map<Double,Double> combineSearchProfiles(
			CompoundParent parentCompound)
			{
		Map<Double,Double> result = new TreeMap<Double,Double>();

		for (Compound species : parentCompound.getCompounds())
		{
			for (Double retentionTime : species.getProfileInSearchWindow().keySet())
			{
				if (result.containsKey(retentionTime))
				{
					Double oldValue = result.get(retentionTime);
					result.put(retentionTime,oldValue
							+ species.getProfileInSearchWindow().get(retentionTime));
				}
				else
				{
					result.put(retentionTime,
							species.getProfileInSearchWindow().get(retentionTime));
				}
			}
		}
		return result;
			}

	public static Set<CompoundParent> computeAbundances(String pathToMzFile,
			Set<CompoundParent> compounds,String pathToSaveFile)
			{
		TargetedAbundanceAnalysisImpl.LOG.info("Computing abundaces for "
				+ pathToMzFile);

		MSXMLParser parser = new MSXMLParser(pathToMzFile);
		Map<Double,Integer> timeScanMap =
			MzXmlUtils.getTimeToScanMap(pathToMzFile,parser);
		Map<Integer,MassScan> retrievedScans = new HashMap<Integer,MassScan>();

		HeatMapBuilder heatmapBuilder =
			new HeatMapBuilder(pathToMzFile,pathToSaveFile);

		double globalMaxMz = Double.MIN_VALUE;
		int globalMinScan = Integer.MAX_VALUE;
		int globalMaxScan = Integer.MIN_VALUE;
		double middleTime = Double.NaN;
		double experimentTotal = 0.0;
		for (CompoundParent compoundParents : compounds)
		{
			// Looking for max in widow
			Map<Scan,Double> maxSearch = new HashMap<Scan,Double>();

			for (Compound species : compoundParents.getCompounds())
			{
				middleTime = species.getRetentionTime();

				List<Integer> scansIndecies =
					MzXmlUtils.getScansFromRetentionTimes(timeScanMap,
							middleTime
							- TargetedAbundanceAnalysisImpl.rtSearchWidth,
							middleTime
							+ TargetedAbundanceAnalysisImpl.rtSearchWidth);

				double minMz =
					species.getMz()
					- TargetedAbundanceAnalysisImpl.mzSumWidthLeft;
				double maxMz =
					species.getMz()
					+ TargetedAbundanceAnalysisImpl.mzSumWidthRight;

				if (maxMz > globalMaxMz)
				{
					globalMaxMz = maxMz;
				}

				for (Integer index : scansIndecies)
				{
					if (!retrievedScans.containsKey(index))
					{
						retrievedScans.put(index,
								new MassScan(parser.rap(index)));
					}
					Scan scan = retrievedScans.get(index).getScan();
					heatmapBuilder.addScan(index,retrievedScans.get(index));
					double intensity =
						MzXmlUtils.getScanSliceIntesitySum(
								retrievedScans.get(index),minMz,maxMz);

					species.addToSearchProfile(
							MzXmlUtils.parseRetentionTime(scan.getRetentionTime()),
							intensity);

					if (intensity > 0.0)
					{
						maxSearch.put(scan,intensity);
					}
				}
			}

			if (maxSearch.size() == 0)
			{
				LOG.warn("Empty Search");
			}

			Scan middleScan =
				TargetedAbundanceAnalysisImpl.getMaxIntensityScan(maxSearch);
			if (middleScan == null)
			{
				LOG.warn("Could not find scan data for molecule " +
						compoundParents.toString());
				
				// Need to 0-out
				for (Compound species : compoundParents.getCompounds())
				{
					species.setAbundance(0.0);
					species.setNumberOfScansInSummationWindow(0);
				}
				
				continue;
			}
			double middleRetentionTime =
				MzXmlUtils.parseRetentionTime(middleScan.getRetentionTime());
			double minRetentionTime =
				middleRetentionTime - TargetedAbundanceAnalysisImpl.rtSumWidth;
			double maxRetentionTime =
				middleRetentionTime + TargetedAbundanceAnalysisImpl.rtSumWidth;

			List<Integer> scansIndecies =
				MzXmlUtils.getScansFromRetentionTimes(timeScanMap,
						minRetentionTime,maxRetentionTime);

			Integer maxScan = Collections.max(scansIndecies);
			Integer minScan = Collections.min(scansIndecies);

			for (Compound species : compoundParents.getCompounds())
			{

				if (minScan <= 0)
				{
					continue;
				}

				if (maxScan > globalMaxScan)
				{
					globalMaxScan = maxScan;
				}

				if (minScan < globalMinScan)
				{
					globalMinScan = minScan;
				}

				heatmapBuilder.addHeatMapTarget(species,minScan,maxScan
						- minScan);
				double sliceTotal = 0;
				species.setNumberOfScansInSummationWindow(maxScan - minScan + 1);
				for (int i = minScan; i <= maxScan; i++)
				{
					if (!retrievedScans.containsKey(i))
					{
						retrievedScans.put(i,new MassScan(parser.rap(i)));
					}

					double minMz =
						species.getMz()
						- TargetedAbundanceAnalysisImpl.mzSumWidthLeft;
					double maxMz =
						species.getMz()
						+ TargetedAbundanceAnalysisImpl.mzSumWidthRight;

					Scan scan = retrievedScans.get(i).getScan();
					double sum =
						MzXmlUtils.getScanSliceIntesitySum(
								retrievedScans.get(i),minMz,maxMz);
					species.addToComputedProfile(
							MzXmlUtils.parseRetentionTime(scan.getRetentionTime()),
							sum);

					sliceTotal += sum;
				}
				double moleculeSum = 0;
				moleculeSum = sliceTotal / species.getP();
				species.setAbundance(moleculeSum);
				experimentTotal += species.getAbundance();
			}
		}
		for (CompoundParent compoundParents : compounds)
		{
			for (Compound species : compoundParents.getCompounds())
			{
				species.setRelativeAbundance(species.getAbundance()
						/ experimentTotal);

				IntensityValue intensityValue = new IntensityValue();
				intensityValue.setNumberOfScans(species.getNumberOfScansInSummationWindow());
				intensityValue.setRawIntensity(species.getAbundance());
				intensityValue.setRelativeIntensity(species.getAbundance()
						/ experimentTotal);
				intensityValue.setFileName(new File(pathToMzFile).getName());

				species.addIntensityValue(intensityValue);
			}
		}

		// Do heatmap
		heatmapBuilder.setMaxMz(globalMaxMz);
		heatmapBuilder.setMinScanToPlot(globalMinScan);
		heatmapBuilder.setMaxScanToPlot(globalMaxScan);
		heatmapBuilder.run();

		return compounds;
			}

	/**
	 * Find the Scan with the highest intensity
	 * @param sums
	 * @return
	 */
	private static Scan getMaxIntensityScan(Map<Scan,Double> sums)
	{
		double max = Double.MIN_VALUE;
		Scan result = null;
		//	System.out.println(sums);

		for (Scan scan : sums.keySet())
		{
			if (sums.get(scan) > max)
			{
				max = sums.get(scan);
				result = scan;
			}
		}
		if (result == null)
		{
			LOG.info("Returning a null for a max scan");
			LOG.info("Number of scans in the map = " + sums.size());
			for (Scan scan : sums.keySet())
			{
				System.out.println(sums.get(scan));
			}

		}

		return result;
	}

	public static void setParams(WidthPair mzSum,WidthPair rtSum,
			WidthPair rtSearch)
	{
		TargetedAbundanceAnalysisImpl.rtSearchWidth = rtSearch.left;
		TargetedAbundanceAnalysisImpl.rtSumWidth = rtSum.left;
		TargetedAbundanceAnalysisImpl.mzSumWidthLeft = mzSum.left;
		TargetedAbundanceAnalysisImpl.mzSumWidthRight = mzSum.right;
	}
}
